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In the "Notes and Discussion" article Ref. [l| the au- 
thors present a very interesting, although purely math- 
ematical, matrix representation of some special polyno- 
mials. Here we apply this mathematical formalism to a 
physical case, namely non-rectilinear propagation of rays 
of light in arbitrary non-homogeneous media. 

Rays of light propagate along rectilinear trajectories 
in air. Therefore, at the generic position x a ray 
(represented by the linear function f{x) = a -\- bx) 
is completely determined by a pair numbers solely: 
f{x) and f'{x). Such a pair may be represented in a 
vector-like form as follows: 



f(x) 



/(^) 
/'(^) 



(1) 



The simple linear relation existing between f(a;i) 
and f(x2) at two arbitrarily chosen positions xi and 
X2 = xi -\- L, with L > 0, is usually written in optics 
textbooks @t3 in the following matrix form: 



/(^2) 



1 L 
1 



/(^i) 



(2) 



The 2x2 matrix in Eq.([2]) is known as the ABCD matrix 
of the optical system (free space, in the present case) and 
fully characterizes the propagation of rays of light in it. 

Literature is rich of examples of ABCD matrices for 
more complicated optical systems as lenses, planes and 
curved dielectric interfaces, mirrors, inhomogeneous me- 
dia with a quadratic index profile et cetera, and combi- 
nations thereof [1, 0] • 

For complex inhomogeneous media, the trajectory of a 
ray of light is more complicated than a straight line and 
cannot be represented anymore by a linear function, ex- 
cept then for a very short distance X2 ~ xi. If this is not 
the case, then one may either subdivide the bulk material 
of macroscopic thickness L in N ^ 1 slices of infinites- 
imal thickness dL = L/N and use N different ABCD 
matrices (one for each slice) [1], or use the alternative 
approach that will be the subject of this note. 



This method is based on the observation that the 
ABCD matrix formalism can be viewed as a lineariza- 
tion of the trajectory of a ray of light around the ini- 
tial point xi [qI, [io| . Such linearization, namely a Tay- 
lor expansion truncated up to and including first order 
terms, is physically meaningful only for X2 close enough 
to xi: X2 ~ xi -I- L/N. But what if X2 is no longer 
close to xi? Does the first order Taylor expansion break 
down? If so, can such expansion be suitably extended? If 
higher order terms must be retained, what is their phys- 
ical meaning? To answer these questions, we put on rig- 
orous basis this linearization procedure showing that a 
2x2 ABCD matrix is a principal sub-matrix [ll| of an 
effective oo x oo matrix describing the full nonlinear dy- 
namics of a curvilinear ray of light. In doing so, we feel 
that re-expressing and generalizing a standard result in 
classical optics (namely, the ABCD matrix formalism) 
in the more familiar language of linear algebra and ba- 
sic calculus would make this topic, usually reserved to 
graduate students in optics, accessible to a more broad 
audience. 

To begin with, let us first re-derive Eq.([2]) for an 
arbitrary linear function f{x) that now we write in the 
following manner: 



f{x) = flo + aix = Oq -I- a\{x — xi), 



(3) 



where Xi is an arbitrary point belonging to the domain of 
the function f{x). If we choose Xi = Q then we retrieve 
the previous expression f{x) — a + bx with ao = a and 
fli = 6. However, for Xi ^ 0, the last equality in Eq.Q 
gives: 



Oo = ao + aiXt = f{xt), 
a\ = ai f'{xi), 



(4a) 
(4b) 



where, for the sake of simplicity, we have introduced the 
notation — ak{xi). 

Since the point Xi is arbitrarly chosen, we can choose 
a different point Xj > Xi and write: 



/(x) = Aq -I- a\{x — Xi) = al + a{{x - Xj). 



(5) 



'Electronic address: ' marco.oniigotti@mpl.mpg.de| 



By equating the factors with the same powers of x at 
the second and third terms in the equation above, we 



2 



obtain Oq = ag + {xj — Xi)a\ and a{ = a\. This can be 
rewritten in the foUowing matrix form; 







1 Xj Xi 










1 







(6) 



This resuh is equivalent with the one written in Eq. ([2|) , 
if we identify ag = f{x2), a{ = f'{x2), flg = f{xi) and 
a\ = f'{xi). The formal derivation of Eq.© via the 
steps (3-6) it is highly redundant for the linear-function 
case. However, it has the virtue to be generalizable to 
the case of non-rectilinear ray propagation. 

Now, in order to describe a ray that propagates in 
an inhomogeneous medium we need a generic smooth 
non-linear function ,f{x) which can be expanded in a 
Taylor series around x — Q 



This expression can be turned into a recursive relation 
by equating terms with the same power of x. Then, for 
k — we have: 



5](-ira>r = E(-i)"«"^"' (12) 

n=0 n=0 

which can be rewritten, after isolating the n = term, as: 



n=l 

For k — \ the same procedure gives 



f{x) = ao H- aix + a2X^ + 



(7) 



For any S M we can write x = x — Xi -\- Xi and insert 
this relation into Eq.© to obtain 



f{x) ^ ao + ai{x - Xi + Xi) + a2{x - Xi -\- Xi)^ -\- ■ ■ ■ 

oo 

= E<(^-^0", (8) 

n=0 

where the coefhcient are given by: 



. 1 a^f{x) 



n\ dx^ 



E 



a-kXi 



(9a) 
(9b) 



Now we can repeat the same reasoning that lead to 
Eq.® and write the following equality: 



E <ix - = E ""(^ " ^j)"' 



(10) 



n=0 



n=0 



where again Xi and Xj are two arbitrarily points on the 
real axis with Xj > Xi. 

By expanding both sides of this equation with the 
help of the Newton's binomial formula one gets 



i = < + E(-i)""'(«"' - <^r')- (14) 



n=2 



This method can be iterated by choosing A: = 2, 3, • ■ • , n 
to generate the following recursive relation: 



«i=«i+ E Q(-i)"-'(«-'-<^r')' (15) 

n=k+l ^ ^ 

with fc = 0, 1, • • • , n. The equation above can be seen as 
a linear algebraic system relating the variables aj^ to the 
quantities a:^. This result can be then written in matrix 
form. Let us introduce the following notation 



/I \n-k n-k 

kr ' ^ 



(16) 



Note that 6'^(fc) = 1. If we now assume to truncate the 
summation in Eq. (jl5p to a finite value N in order to be 
able to represent these results (formally described by 
an oo X oo matrix) with a finite dimensional matrix as 
follows: 



N 

E 

n—k 



N 

bi{k)ai^Y.^n{k)< 

n—k 



(17) 



E < E 



n=0 



k=0 



\n — k 



n=0 k=0 



E<E : -'(--.)"-'=■ (11) 



By then defining — (oq, aj, • ■ • , a]^) and 

a' = (oq, a^, • • • , o^y), it is possible to rewrite the 
previous equation in matrix form as follows: 
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61(0) biiO) ...1 








&i(0) 6^(0) 






6i(l) mi) ••• 








h\{l) h\(l) ••• 




a\ 


b{{2) ••• 








61(2) ■•• 







(18) 



If we call B(j) the matrix on the left-hand-side of the 
previous equation and B(i) the one on the right-hand- 
side, the previous relation can be written in the following 
compact form: 

B(j)a^' = B(i)a\ (19) 

where again the shorthand notations B(i) = B(a;i) and 
B(j) = B(a;j) are used for the sake of clarity. Solving 
for by multiplying on the left both sides of the previ- 
ous equation by B(j)~^ and defining A = B~^(j)B(i), 
we can write the relation between the vectors a? and a* as 



This equation has a straightforward physical meaning: 
it illustrates the fact that propagation across two 
consecutive distances Li and L2 can be described as a 
single propagation along the distance Li -I- L2. From a 
mathematical point of view, this is a signature of the 
semigroup property of our generalized ABCD matrices. 
With the help of a suitable mathematical software for 
algebraic manipulation, it is not difficult to verify via 
explicit N X N matrix multiplications, that Ea. (P5|) is 
valid for arbitrary N. Thus, by iteration, one can easily 
convince oneself that the matrix A satisfies the following 
relation !3|: 



Aa' 



(20) 



The matrix A is our sought generalized ABCD matrix, 
whose expression is the following: 



N 



N 



l[A{Ln)=A(Y,Ln 



(24) 



A - B-i(j)B(z) 

rl L L3 .... 

1 2L ••• 

^ 1 3L ••• 

1 ••• 

EE A(i), (21) 

where L = Xj — Xi. Note that this matrix contains the 
usual (i.e. linear) ABCD matrix defined in Eq.(l2]) as 
the first 2x2 principal sub matrix . This sub matrix 
verifies the following equality: 



"1 Ll' 




■ 1 


L2 ' 




' 1 L1+L2' 


1 







1 




1 



(22) 



namely A{Li)A{L2) = A{Li + L2). The next step 
is to extract the principal 3x3 submatrix from Eq. 
(PT|). i.e. to consider the first nonlinear order of the 
approximation of the function f{x). In this case, the 
relation A{Li)A{L2) = A{Li + L2) still holds, as can be 
checked by direct calculation: 



1 Ll Ll 




1 L2 Ll 


1 2Li 




1 2L2 


1 




1 




" 1 Li+L 


2 


{Li+L2f' 




1 


2(Li+i2) 







1 



This result has a straightforward interpretation: the 
propagation of the function through the total distance 
Li + L2 + ■ ■ ■ + Lf^ can be achieved by consecutive prop- 
agation across the distances Li,L2,- • • ,Lm- 

This analogy is not accidental. A closer inspection 
to Eq. ((20t reveals in fact that this equation tells us 
how the value of the function f{x) in a point Xj can be 
calculated knowing the value of the same function in a 
point Xi < Xj. With this in mind, we can calculate the 
derivative of f{x) as follows: 



df{x) fix + Ax) - fix) 

— i — = A 

ax Ax-^o Ax 




= limf^^ — -)a' = -Da\ (25) 

L^O \ L J 

where we choose Ax = Xj — Xi = L so that we can 
express /(x -|- Ax) as a^ and /(x) as a*. Note that this 
does not cause any loss of generality, since the definition 
of derivative involves only the concept of neighboring 
points and, as discussed previously, the quantities a* 
and a^ represent the value of the function /(x) in two 
arbitrary neighboring points. Note moreover that in the 
last equality we used Eg. d^D)) to write a-' as a function 
of a*. Here, D is the matrix representation of the 
differential operator d/dx [l| 
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D 
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1 
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••• 











CO 











••• 



(26) 



It is not difficult to see, via an expHcit calculation, 
that the generalized ABCD matrix is related to the 
differential operator by the following formula: 



OO 

fe=0 



(LD) 



(27) 



From didactic point of view it is very gratifying to 
see that we were able to reproduce the already known re- 



sult for the translation operator acting upon an arbitrary 
function. In fact, Eq. (P7)) gives an actual physical repre- 
sentation of the well-known translation operator e^i^.^'^^^ 
1131, such that: 



/(a;)|x=o 



E 

fc=0 



fc! 



= .f{L). 



(28) 



In our case, in fact, the action of the A matrix com- 
pletely defines the vector F (whose elements are the all 
a\) at the point X2 knowing the expression of the F 
vector at the point xi, i.e. 



F(a;2) 



F(a;i) = A(L)F(xi). 



(29) 
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